%% mainScript.m
% This script computes the baseline model solution and all moments listed
% in Table 2. It also creates Figures 1 and 2 (at least approximately; the
% figures in the paper have been hand-edited).
%% script settings
% the following parameters control the behavior of the solution algorithm
% and display settings

printMoments = true; % if true, prints moments for table 2
createPlots = true; % if true, creates the plots for Figures 1 and 2
% algorithm parameters
%   - pde algorithm
nGridpoints = 800; % number of points for \tilde{\sigma} grid
dt = 0.33; % time step size in pde solution algorithm
nSolutionSteps = 1500; % maximum number of iterations (time steps) in pde solver
convCrit = 5e-6; % convergence criterion for pde solver (based on maximum relative change in function values)
%   - simulation
nSimPeriods = 1000; % number of periods (years) to simulate
nSimBurnIn = 50; % this many periods are discarded when computing moments
simTimeStep = 0.01; % time step in simulation
% plotting parameters
plotsLbSigmaI = 0.15; % lowest value for \tilde{\sigma} in plots (Figures 1 and 2)
plotsUbSigmaI = 0.95; % highest value for \tilde{\sigma} in plots (Figures 1 and 2)

%% fixed parameters
% these are parameters fixed prior to the moment matching exercise

% sigma process is externally calibrated
sigma = sqrt(0.29);  % mean reversion level of \tilde{\sigma}
stateVol = 0.4; % volatility in state evolution (state is \tilde{\sigma}^2)
stateMeanReversion = 0.67;% % mean reversion speed in state evolution (state is \tilde{\sigma}^2)

% skin-in-the-game constraint is externally calibrated
chi = 0.3;

% depreciation rate does not really matter/is exeternally calibrated
delta = 0.055;

%% variable parameters
% these have been varied manually to achieve the calibration targets
% (approximately)

% direct parameters
gamma = 6;
phi = 8.1;

% linear rule parameters
alphaA_rel = 0.115;
alphaBondGrowth = 0.12;

% target parameters
% (Transformed variants of model parameters that are more closely linked to
% certain target moments. The link is precise in steady state conditional
% on matching \vartheta but only approximate otherwise. The function 
% preCalibration transforms these variables into model parameters.)
debtGdp = 0.7;
capitalGdp = 3.5;
investmentGdp = 0.2;
govSpendingGdp = 0.22;
surplusGdp = -0.0018;
invRate = 0.125;

%% compute model solution, moments, baseline simulation
% pack all algorithm parameters in a struct array
algParams = aux.pack(nGridpoints,dt,nSolutionSteps,convCrit,nSimPeriods,nSimBurnIn,simTimeStep);

% determine remaing model parameters using pre-calibration
[a0,g0,rho,adjBondGrowth0,iota0] = preCalibration(debtGdp,capitalGdp,investmentGdp,govSpendingGdp,surplusGdp,invRate,phi);

% compute absolute sensitivities in linear rules (inputs are relative to
% level except for policy whose level is close to zero)
alphaA = alphaA_rel*a0;

% pack all model parameters in a struct array
params = aux.pack(sigma,stateVol,stateMeanReversion,a0,g0,iota0,rho,chi,gamma,phi,...
  adjBondGrowth0,delta,alphaA,alphaBondGrowth);

% call main solution function
[moments,modelSolution,modelSim,xGrid,xDrift,xVolatility] = computeMomentsFromParams(params,algParams);

%% plot for equilibrium dynamics (Figure 1)

plotsUseIdx = modelSolution.sigmaI >= plotsLbSigmaI & modelSolution.sigmaI <= plotsUbSigmaI;
plotsGridSigmaI = modelSolution.sigmaI(plotsUseIdx);

% stationary density

stateMeanHeston = sigma^2;
gammaParamAHeston = 2*stateMeanReversion*stateMeanHeston/stateVol^2;
gammaParamBHeston = stateVol^2/2/stateMeanReversion;
quantilesHeston = icdf('Gamma',[0.0005,0.9995],gammaParamAHeston,gammaParamBHeston);

tempGrid = linspace(quantilesHeston(1),quantilesHeston(2),10000);
tempSigmaI = sqrt(tempGrid);
tempTransfFactor = 2*sqrt(tempGrid); %factor from integral transformation formula
tempStatDensity = pdf('Gamma',tempGrid,gammaParamAHeston,gammaParamBHeston).*tempTransfFactor;

densityScaling = 0.98*(max(modelSolution.qB(plotsUseIdx)) - 0.93*min(modelSolution.qB(plotsUseIdx)))/max(tempStatDensity);
densityOffset = 0.93*min(modelSolution.qB(plotsUseIdx));

figure(1);
yyaxis left;
area(tempSigmaI,densityOffset + densityScaling*tempStatDensity,'FaceColor','#cccccc','EdgeColor','#808080');
hold on
plot(plotsGridSigmaI,modelSolution.qB(plotsUseIdx),'LineStyle','-','LineWidth',2);
xlim([plotsGridSigmaI(1),plotsGridSigmaI(end)]);
ylim([floor(min(modelSolution.qB(plotsUseIdx))*10)/10,ceil(max(modelSolution.qB(plotsUseIdx))*10)/10]);
ylabel('$q^B$','Interpreter','latex');
yyaxis right;
plot(plotsGridSigmaI,modelSolution.qK(plotsUseIdx),'LineStyle','-','LineWidth',2);
ylim([floor(min(modelSolution.qK(plotsUseIdx))*10)/10,ceil(max(modelSolution.qK(plotsUseIdx))*10)/10]);
ylabel('$q^K$','Interpreter','latex');

xlabel('$\tilde{\sigma}$','Interpreter','latex');
ax = gca; %axes handle
ax.FontSize = 14;
ax.TickLabelInterpreter = 'latex';
ax.Box = 'on';

%% compute additional pde solutions
% the following function solves two additional sets of pdes:
%  - pdes for bond value decomposition (cash flow and service flow terms)
%    in the baseline model
%  - pdes for model without bonds (for flight-to-safety volatility)

[pvSurplusesConsNum,pvRiskSharingConsNum,logValueNoBonds,sigmaTotalWealthNoBonds,aggrPriceOfRiskNoBonds] ...
    = solveModelAdditionalPdes(xGrid,xDrift,xVolatility,modelSolution,params,algParams);
  
% numbers for flight-to-safety volatility (stated in Section 4.5

differentiator = finiteDifferentiator(xGrid,3);
amg = modelSolution.a-modelSolution.govSpending; %productivity (a) minus (m) government spending (g)
sigmaAmg = differentiator.d1(log(amg)).*xVolatility;
noDebtSigmaQ = phi*amg./(1+phi*amg).*sigmaAmg;
excessReturnVolatilityNoDebt = -noDebtSigmaQ;
excessReturnVolatilityNoDebtInterpolant = griddedInterpolant(xGrid,excessReturnVolatilityNoDebt);
simNoDebtExcessReturnVolatility = excessReturnVolatilityNoDebtInterpolant(modelSim.x);
meanExcessReturnVolatilityNoDebt = mean(simNoDebtExcessReturnVolatility);
meanExcessReturnVolatilityBaseline = moments.meanEquityPremium/moments.sharpeRatio;


%% plot decomposition figure (Figure 2)

figure(2);
hold on
plot(plotsGridSigmaI,[pvSurplusesConsNum(plotsUseIdx),pvRiskSharingConsNum(plotsUseIdx)],'LineWidth',2);
plot(plotsGridSigmaI,modelSolution.qB(plotsUseIdx),'k--','LineWidth',2);
xlim([plotsGridSigmaI(1),plotsGridSigmaI(end)]);
ylabel('$q^{B,CF}, q^{B,SF}, q^B$','Interpreter','latex');
xlabel('$\tilde{\sigma}$','Interpreter','latex');
legend({'$q^{B,CF}$', '$q^{B,SF}$', '$q^B$'},'Interpreter','latex');
plot(plotsGridSigmaI,zeros(size(plotsGridSigmaI)),'k:');
ax = gca; %axes handle
ax.FontSize = 14;
ax.TickLabelInterpreter = 'latex';
ax.Box = 'on';

%% print moments
if ~printMoments
  return
end

fprintf('-------------------------------------------------\n')
fprintf('######    moments from model simulation:   ######\n')
fprintf('-------------------------------------------------\n')
fprintf('  targeted moments\n')
fprintf('-------------------------------------------------\n')
fprintf('std(Y):          %.4f  -- target:  %.4f\n',moments.stdOutput,0.013)
fprintf('std(C)/std(Y):   %.4f  -- target:  %.4f\n',moments.stdCons/moments.stdOutput,0.64)
fprintf('std(I)/std(Y):   %.4f  -- target:  %.4f\n',moments.stdInv/moments.stdOutput,3.38)
fprintf('std(S/Y):        %.4f  -- target:  %.4f\n',moments.stdSurplusOutputRatio,0.011) 
fprintf('E[C/Y]:          %.4f  -- target:  %.4f\n',moments.meanConsOutpRatio,0.56)
fprintf('E[G/Y]:          %.4f  -- target:  %.4f\n',moments.meanGovSpendingOutpRatio,0.22)
fprintf('E[S/Y]:         %.4f  -- target: %.4f\n',  moments.meanSurplusOutputRatio,-0.0005) 
fprintf('E[I/K]:          %.4f  -- target:  %.4f\n',moments.meanIota,0.12)
fprintf('E[q^K K/Y]:      %.4f  -- target:  %.4f\n',moments.meanCapitalOutpRatio,3.73)
fprintf('E[q^B K/Y]:      %.4f  -- target:  %.4f\n',moments.meanDebtOutpRatio,0.71)
fprintf('E[r^e-r^b]:      %.4f  -- target:  %.4f\n',moments.meanEquityPremium,0.034)
fprintf('sharpe ratio:    %.4f  -- target:  %.4f\n',moments.sharpeRatio,0.31);
fprintf('-------------------------------------------------\n')
fprintf('  untargeted moments\n')
fprintf('-------------------------------------------------\n')
fprintf('corr(C,Y):       %.4f\n',moments.corrConsOutput)
fprintf('corr(I,Y):       %.4f\n',moments.corrInvOutput)
fprintf('corr(S/Y,Y):     %.4f\n',moments.corrSurplusOutputRatioOutput)  
fprintf('std(q^B K/Y):    %.4f\n',moments.stdDebtOutputRatio)  
fprintf('E[r^f]:          %.4f\n',moments.meanRiskFreeRate)  
fprintf('std(r^f):        %.4f\n',moments.stdRiskFreeRate)  
fprintf('-------------------------------------------------\n')
fprintf('  other (numbers stated in Section 4.5)\n')
fprintf('-------------------------------------------------\n')
fprintf('std(r^e-r^b) with bonds:    %.4f \n',meanExcessReturnVolatilityBaseline)
fprintf('std(r^e-r^b) without bonds: %.4f \n',meanExcessReturnVolatilityNoDebt);
fprintf('-------------------------------------------------\n')
